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Abstract. We present a lattice Boltzmann model to describe the spread- 
ing of droplets on topologically patterned substrates. We apply it to 
model superhydrophobic behaviour on surfaces covered by an array of 
micron-scale posts. We find that the patterning results in a substantial 
increase in contact angle, from 110° to 156°. 

1 Introduction 

A droplet in contact with a substrate will try to spread to an equilibrium shape 
determined by Young's law which describes the balance of surface tensions. There 
are many parameters which affect this process. For example surface disorder in 
the form of chemical or topological heterogeneities can pin a droplet or change 
its final shape. This has usually been viewed as a nuisance in experiments and 
applications. However with the advent of microfabrication techniques it is be- 
coming possible to harness controlled surface topologies to explore new physical 
phenomena. 

A beautiful example of this, inspired by the leaves of the lotus plant, is a 
superhydrophobic substrate. The angle 9 between the tangent plane and the 
droplet is usually known as the contact angle. The higher the contact angle the 
more repellent the surface. There are applications, for example raincoats and 
windscreens, where repellent surfaces are highly desirable. Surface coatings and 
chemical modifications of the substrate are common ways to increase the contact 
angle but it is difhcult to achieve an angle of more than 120°. However surfaces 
patterned with posts on a micron length scale allow contact angles of 160° to be 
reached [1,2]. 

Another example is nanopatterning of substrates which is of particular in- 
terest for the semiconductor industry. Applying a microstamp on a thin heated 
polymer film coating a substrate can create nanoscale patterns. The film dewets 
along the faces of the mold to form either thin lines or small patches. The present 
model has been used to simulate this new technique [3]. 

The aim of this paper is to present a lattice Boltzmann algorithm which 
can be used to investigate the behaviour of droplets on topologically patterned 
substrates. Lattice Boltzmann is a particularly appropriate approach in that it 
solves the Navier Stokes equations but also inputs the thermodynamic infor- 
mation such as surface tensions needed to describe the behaviour of droplets. 



Moreover its natural length scale, for fluids such as water, is of order microns 
where much of the exciting new physics is expected to appear. The method has 
already shown its capability in dealing with spreading on surfaces with chemical 
patterning [4]. 

In section 2 we summarise the algorithm and, particularly, describe the new 
thermodynamic and velocity boundary conditions needed to treat surfaces with 
topological patterning. In section 3 we present results for a substrate patterned 
by an array of posts. The patterning leads to a considerable increase in contact 
angle. Finally we discuss directions for future work using this approach. 

2 The mesoscopic model 

We consider a liquid-gas system of density n(r) and volume V. The surface of 
the substrate is denoted by S. The equilibrium properties are described by the 
free energy 

ipb{n) is the free energy in the bulk. We choose a Van der Waals form 

Mn) = Pc {^n + 1)' {^1 - 2i/„ + 3 - 2/3r^) (2) 

where = (n — n^jnc^ t^, = (T,, — T)/Tc and pc = 1/8, ric = 7/2 and Tc = 4/7 
are the critical pressure, density and temperature respectively and /? is a constant 
typically equal to 0.1. The bulk pressure 

Pb = Pc{vn + ifiii^l - 2i/„ + 1 - 2/3r™). (3) 

The derivative term in equation (1) models the free energy associated with 

an interface, k, is related to the surface tension, ipdng) = — 01*^8 H is the 

Cahn surface free energy [5] which controls the wetting properties of the fluid. 

The lattice Boltzmann algorithm solves the Navier-Stokes equations for this 
system. Because interfaces appear naturally within the model it is particularly 
well suited to looking at the behaviour of moving drops. 

2.1 The Lattice Boltzmann algorithm 

The lattice Boltzmann approach follows the evolution of partial distribution 
functions fi on a regular, rf-dimensional lattice formed of sites r. The label i 

denotes velocity directions and runs between and z. DdQz + 1 is a standard 
lattice topology classification. The D3Q15 lattice we use here has the follow- 
ing velocity vectors Vi-. (0,0,0), (±1,±1,±1), (±1,0,0), (0,±1,0), (0,0, ±1) in 
lattice units as shown in fig. 1. 

The lattice Boltzmann dynamics are given by 



Mr + Atvu t + At) = Mr, t) + - (//« (r, t) - /,(r, t)) 



(4) 
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Fig. 1. Topology of a D3Q15 lattice. The directions i are numbered and correspond to 
the velocity vectors Vj. 



where At is the time step of the simulation, r the relaxation time and f^'^ the 
equilibrium distribution function which is a function of the density n = J2i=o fi 
and the fluid velocity u, defined through the relation 



nu 

i=0 



The relaxation time tunes the kinematic viscosity as [6] 

Ar^ d , 1 



V = 



At 2 



(t-t.) (6) 



where Ay is the lattice spacing and C2 and C4 are coefficients related to the 
topology of the lattice. These are equal to 3 and 1 respectively when one considers 
a DSQIS lattice (see [7] for more details). 

It can be shown [8] that equation (4) reproduces the Navier-Stokes equations 
of a non-ideal gas if the local equilibrium functions are chosen as 

fr =^CT + B^UaVia + CaM^ + D^UaUffViaVi^ + Gaa^ViaVip, i > 0, 

/o'="-E/r (7) 
i=l 

where Einstein notation is understood for the Cartesian labels a and /3 (i.e. 
ViaUa = X^c« ViaUa) and where a labels velocities of different magnitude. A 
possible choice of the coefficients is [9] 



_ WcrTl _ WgU _ SWaU 



2c2 ' " 2c4 



1 



^177 ~ ('^(^7^)^ + ^UU-yd-yTl) , G277 = 0, 



G^-yS = {K{d^n){dsn) + v{u^dsn + ugdju)) (8) 
where wi = 1/3, W2 = 1/24 and c = Ar/At. 



2.2 Wetting boundary conditions 

The major challenge in dealing with patterned substrates is to handle the bound- 
ary conditions correctly. We consider first wetting boundary conditions which 
control the value of the density derivative and hence the contact angle. For flat 
substrates a boundary condition can be established by minimising the free energy 

(1) [5] 

s.Vn = -— (9) 

K 

where s is the unit vector normal to the substrate. It is possible to obtain an 
expression relating to the contact angle 9 as [9] 



= 2l3Tyj^2pcK sign (| - ^) y cos ^ (l - cos ^) (10) 

where a = cos~^(sin^ 9) and the function sign returns the sign of its argument. 

Equation (9) is used to constrain the density derivative for sites on a flat part 
of the substrate. However, no such exact results are available for sites at edges 
or corners. We work on the principle that the wetting angle at such sites should 
be constrained as little as possible so that, in the limit of an increasingly fine 
mesh, it is determined by the contact angle of the neighbouring flat surfaces. 

For edges (labels 9 — 12 in fig. 2) and corners (labels 1 — 4) at the top of 
the post each site has 6 neighbours on the computational mesh. Therefore these 
sites can be treated as bulk sites. 

At bottom edges where the post abuts the surface (labels 13 — 16 in fig. 2) 
density derivatives in the two directions normal to the surface (e.g. x and z for 
sites labeled 13) are calculated using 

d,n = d^/yn=-^^ (11) 

where the middle term constrains the density derivative in the appropriate di- 
rection X or y. 

At bottom corners where the post joins the surface (labels 5 — 8 in fig. 2) 
density derivatives in both the x and y directions are known. Therefore these 
sites are treated as planar sites. 



2.3 Velocity boundary conditions 



We impose a no-slip boundary condition on the velocity. Because the collision 
operator (the right hand side of equation (4)) is applied at the boundary the 
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Fig. 2. Sketch of a post on a substrate. Encircled numbers label sites in different 
topological positions. Labels 26 and 27 denote sites on the bottom {z = Zmin) and the 
top {z = Zmax) of the domain respectively. 

usual bounce-back condition is not appropriate as it would not ensure mass 
conservation [10]. 

Indeed after applying equation (4) there are missing fields on the substrate 
sites because no fluid has been propagated from the solid. Missing fields are 
determined to fulfill the no-slip condition given by equation (5) with u = 0. This 
does not uniquely determine the /i's. For most of the cases (i.e. 1 — 20) arbitrary 
choices guided by symmetry are used to close the system. This is no longer 
possible for sites 21 — 27 where four asymmetrical choices are available. Selecting 
one of those solutions or using a simple algorithm which chooses one of them at 
random each time step leads to very comparable and symmetrical results. Hence 
we argue that an asymmetrical choice can be used. Possible conditions, which 
are used in the results reported here, are listed in table 1. 

The conservation of mass is ensured by setting a suitable rest field, /o, equal 
to the difference between the density of the missing fields and the one of the 
fields entering the solid after collision. 

3 Results 

As an example we consider here the superhydrophobic behaviour of droplet 
spreading on a substrate patterned by square posts arranged as in fig. 3. 

The size of the domain = 80 X 80 X 80 and the height, spacing 

and width of posts are h = 5, d = 8 and w = 4 respectively. A spherical droplet 
of radius i? = 30 is initially centered around the point {x; y; z) = (41; 41; 36). The 
contact angle Oinput = 110° is set on every substrate site. The surface tension 
and the viscosity are tuned by choosing parameters k = 0.002 and r = 0.8 



Label 


Conditions 


1 




2 


fr =/8 


3 


/g = /lo 


4 


/ll = /l2 


5 


/S =/6 

fl -i — ( h — fd — fl + fr>)/2+ fa 

^13 W3 J4 Jl—J2)j^^J\j 

+ /l4 — /lO 
/ll = (/l - /2)/2 - /9 + /lO + /l2 
/7 = (-/3 + /4)/2 + /s - /g + /lO 


6 


/S =/6 

fl, = f fo — f/y2 — fl1 + fl9 + fl4 
/9 = (/l - /2)/2 - /ll + /lO + /l2 
fr = {-fl + /2 - /s + /4)/2 + /8 
+/ll - /l2 


7 


/5 =/6 

/ll = (/3 - /4)/2 - /l3 + /l2 + /l4 
/) = {fl - /2 - /a + /4)/2 + /l3 
^/l4 + /lO 

fr = {-fl + /2)/2 + /8 - /la + /i4 


8 


/5 =/6 

/ll = (/a - /4 + /i - /2)/2 + fr 

-fs + /l2 
/9 = (-/3 + /4)/2 - /7 + /8 + /lO 
/l3 =(-/l+/2)/2-/V+/8+/l4 


9 


/l3 = /l4 


/7 = /8 


10 


/g = /lo 


/7 = /s 


11 


/s = /i() 


/ll = /l2 


12 


/l3 = /l4 


/ll = /l2 


13 


/5 =/6 

/l =2(-/io+/9+/ll-/l2)+/2 
/l3 = (/3 - /4)/2 - /ll + /l2 + /l4 
fr = (-/3 + /4)/2 + /8 - /g + /lO 


14 


/5 =/6 

/g = (/i - /2)/2 - /ll + /lo + /i2 

fr =(-/l+/2)/2 + /8-/l3+/l4 
/3 = 2(-/l2 + /ll + /l3 - /l4) + /4 


15 


/5 =/6 

/2 = 2(-/l4 + /7 + /l3 - /8) + fl 
/ll = (/3 - /4)/2 - /l3 + /l2 + /l4 

/g = (-/a + /4)/2 + /8 - A + /lo 



Label 


Conditions 


16 


fr, = /(i 

/ll =(/l-/2)/2-/9+/l0+/l2 
/l3 = {-fl + /2)/2 - /7 + /8 + /l4 
/4 = 2(-/io + /7 + /9 - /s) + /a 


17 


f I (1 = / '1 ! / I = / "1 4 


18 


/7 =/« : /i2 = /n 


19 


/g = /lo ; /l4 = /l3 


20 


f» = fr ; /ll = /i2 


21 


fl ~ fi : fr — fs 
fi2 = (-/a + /4)/2 + /ll 

/l3 = (-/5 + /6)/2 + /l4 

/l« = (/3 - /4 + /s - /6)/2 + /g 


22 


JO — Ta ' 7"^ — To 

^3 J4 , J 7 — J 8 

/9 = (-/S + /6)/2 + /lO 

/l4 = {fl - /2 + /5 - /6)/2 + /l3 

/l2 = {-fl + /2)/2 + /ll 


23 


f2 = fl ; /8 = /7 

/ll = (/3 - /4)/2 + /l2 
/l4 = (/6-/6)/2 + /l3 

/9 =/l() + (-/3 + /4-/5+/6)/2 


24 


Ta • jA To 
7" J' 1 J ^ J 

/ll =(/l-/2)/2 + /l2 

/lO = (/s - /6)/2 + /g 

/l3 = (-/S + /6 - /l + /2)/2 + /l4 


25 


fr =/8 ; /5=/6 

/ll = (/l - /2 + /3 - /4)/2 + /12 

/9 = (-/3 + /4)/2 + /lO 

/l3 =(-/l+/2)/2 + /l4 


26 


fr = fs ; /s = /e 

/ll = (/l - /2 + /3 - /4)/2 + /l2 

/. = (-/3 + /4)/2 + /l„ 
/i:'. = {-fl T ,/2)/2 ^ ,/i4 


27 


/e = /s ; fs ^ fr 

fio = (/3 - /4)/2 + /g 

/l4 = {fl - /2)/2 + /l3 

/l2 = -(-/3 + /4 - /l + /2)/2 + /ll 



Table 1. Velocity boundary conditions. 
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Fig. 3. Sketch of the substrate. Dimensions are expressed in simulation units. Shaded 
areas are posts. 



respectively. The liquid density ni and gas density Ug are set to n; = 4.128 and 
Ug = 2.913 and the temperature T = 0.4. 

Fig. 4 shows the final state attained by the droplet for different substrates 
and initial conditions. For comparison fig. 4(a) shows a planar substrate. The 
equilibrium contact angle is 9a = 110° = dmput as expected [9]. In fig. 4(b) the 
substrate is patterned and the initial velocity of the drop is zero. Now the contact 
angle is 6b = 156°, a demonstration of superhydrophobic behaviour. Fig. 4(c) 
reports an identical geometry but a drop with an initial impact velocity. Now the 
drop is able to collapse onto the substrate and the final angle is 0b = 130°. These 
angles are compatible with the ones reported in [2] where similar parameters are 
considered. 

For the parameter values used in these simulations the state with the droplet 
suspended on the posts has a slightly higher free energy than the collapsed state. 
It is a metastablc state and the droplet needs an impact velocity to reach the 
true thermodynamic ground state. For macroscopic drops gravity will also be 
important in determining whether the drop remains suspended on top of the 
posts. Extrand has predicted the minimum post perimeter density necessary 
for a droplet to be suspended [11]. A next step will be to add gravity to the 
simulation to compare to his prediction. 

Superhydrophobicity occurs over a wide range of d, the distance between the 
posts. For suspended drops of this size and d > 12 the drop resides on a single 
post and the contact angle is 170°. For < 12 the contact angle lies between 




Fig. 4. Final states of a spreading droplet. The right column reports cuts at y = 41. 
(a) The substrate is flat and homogeneous, (b) The substrate is decorated with posts 
and the initial velocity of the droplet is 0. (c) Same geometry as (b) but the droplet 
reaches the substrate with a velocity O.OlAr/ At. Each of these simulations ran for 
approximately 8 hours on 8 processors on a PC cluster. 



148° and 156" with the range primarly due to the commensurabihty between 
drop radius and post spacing. 

It is of course also of interest to look further at the dynamics of the spreading. 
The droplet random motion reported in [2] and the bouncing back of droplets 
on nanotubes [12] pose many interesting directions for future research. 

4 Conclusion 

We have proposed a lattice Boltzmann model describing the spreading of a 
droplet on topologically patterned substrates. As an example, we have consid- 
ered a substrate patterned with posts and found superhydrophobic behaviour in 
agreement with experiments [1,2]. 

The algorithm gives us the capability to explore a wide variety of interesting 
problems concerning droplet behaviour on novel substrates and in microfluidic 
devices. 
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